#ifndef MINCOSTFLOW_H_
#define MINCOSTFLOW_H_

#include <assert.h>
#include <vector>
#include <utility>

#include <boost/noncopyable.hpp>

#include "graph_support/Graph.h"

#include "graph_support/utils/GraphTemplateUtils.h"


namespace GRAPH_SUPPORT_NS
{
namespace alg
{
namespace flow
{

template< typename GraphT >
class MinCostFlow : private boost::noncopyable
{
private:
	typedef GraphT Graph;
	typedef typename Graph::GrVertex GrVertex;
	typedef typename Graph::GrEdge GrEdge;
	typedef typename Graph::GrEdgePtr GrEdgePtr;
	typedef typename Graph::VertexIterator VertexIterator;
	typedef typename Graph::EdgeIterator EdgeIterator;
	typedef typename Graph::EdgeAttributes EdgeAttributes;
	typedef typename EdgeAttributes::Capacity Capacity;
	typedef typename EdgeAttributes::Cost Cost;

public:
	MinCostFlow( Graph & gr )
		: gr( gr )
		, totalFlow( 0 )
		, totalCost( 0 )
	{
	}

	void operator() ( const GrVertex & srcVert, const GrVertex & destVert, const Capacity & requiredFlow )
	{
		eraseFlow();
		runMinCostFlow( srcVert, destVert, requiredFlow );
	}

	Capacity getTotalFlow() const
	{
		return totalFlow;
	}

	Capacity getFlow( const GrVertex & a, const GrVertex & b ) const
	{
		return utils::getFlow( gr, a, b );
	}

	Cost getTotalCost() const
	{
		return totalCost;
	}

	Cost getCost( const GrVertex & a, const GrVertex & b ) const
	{
		std::vector< GrEdgePtr > buf;
		findEdges( gr, a, b, buf );
		Capacity result = std::accumulate( buf.begin(), buf.end(), ( Capacity ) 0,
				boost::bind( std::plus< Capacity >(), _1,
					boost::bind( std::multiplies< Capacity >(),
						boost::bind( & EdgeAttributes::getFlow, boost::bind( & GrEdge::getAttributes,
								boost::bind( & utils::ptr2ref< GrEdge, GrEdgePtr >, _2 ) ) ),
						boost::bind( & EdgeAttributes::getCost, boost::bind( & GrEdge::getAttributes,
								boost::bind( & utils::ptr2ref< GrEdge, GrEdgePtr >, _2 ) ) )
						)
				) );
		return result;
	}

private:
	void eraseFlow()
	{
		// TODO: use std::for_each after ptr2ref - fix
		VertexIterator stopvit = gr.getVerticesEnd();
		for ( VertexIterator vit = gr.getVerticesBeg(); vit != stopvit; ++vit )
		{
			GrVertex vertex = * vit;
			EdgeIterator stopeit = gr.getEdgsEnd( vertex );
			for ( EdgeIterator eit = gr.getEdgsBeg( vertex ); eit != stopeit; ++eit )
			{
				( * eit )->getAttributes().eraseFlow();
			}
		}
		totalFlow = 0;
		totalCost = 0;
	}

	void runMinCostFlow( const GrVertex & srcVert, const GrVertex & destVert, const Capacity & requiredFlow )
	{
		size_t NOT_A_VERTEX_ID = GRAPH_SUPPORT_NS::Graph::NOT_A_VERTEX_ID;

		size_t srcVid = gr.vert2num( srcVert );
		size_t destVid = gr.vert2num( destVert );

		int vertexCount = gr.getVertexCount();
		std::vector< Cost > distBuf( vertexCount );
		std::vector< Cost > fiBuf( vertexCount );
		std::vector< GrEdgePtr > edgeFromParentBuf( vertexCount );
		std::vector< int > visitedBuf( vertexCount, 0 );
		int curVisitedTrue = 1;
		int * visited = & visitedBuf[0];
		std::vector< char > dvisitedBuf( vertexCount, ( char ) false );
		bool * dvisited = ( bool * ) & dvisitedBuf[0];

		bool allNegCyclesCancelled = false;

		// TODO: explicitly neg cyc canc --> potentials
		while ( totalFlow < requiredFlow )
		{
			std::pair< bool, size_t > fordBellmanResult( false, NOT_A_VERTEX_ID );
			if ( allNegCyclesCancelled )
			{
				// TODO: swap fi, dist
				for ( int i = 0; i < vertexCount; i++ )
				{
					fiBuf[i] = distBuf[i];
				}
				fordBellmanResult.first =
						runPotentialDejkstra( srcVid, destVid, & distBuf[0], & fiBuf[0], dvisited, & edgeFromParentBuf[0] );
				for ( int i = 0; i < vertexCount; i++ )
				{
					distBuf[i] = distBuf[i] - fiBuf[srcVid] + fiBuf[i];
				}
				assert( distBuf[srcVid] == 0 );
			} else
			{
				fordBellmanResult =
						runFordBellman( srcVid, destVid, & distBuf[0], & edgeFromParentBuf[0] );
			}
			bool negativeCycleExists = fordBellmanResult.second != NOT_A_VERTEX_ID;
			// TODO: fordBellmanResult.first is not semantically
			if ( ! fordBellmanResult.first && ! negativeCycleExists )
			{
				break;
			}

			size_t startVid = srcVid;
			size_t finishVid = destVid;
			if ( negativeCycleExists )
			{
				assert( ! allNegCyclesCancelled );
				size_t negativeCycleChild = fordBellmanResult.second;
				size_t negativeCyclePart = discoverNegativeCyclePart( negativeCycleChild, & edgeFromParentBuf[0], visited, curVisitedTrue );
				startVid = negativeCyclePart;
				finishVid = negativeCyclePart;
			} else
			{
				allNegCyclesCancelled = true;
			}
			pushViaPath( startVid, finishVid, & edgeFromParentBuf[0], requiredFlow - totalFlow );
		}
	}

	/**
	 * std::pair<
	 * 		bool: if path 4 push found,
	 * 		size_t: if negative cycle exists => any vertex of cycle and NOT_A_VERTEX_ID otherwise
	 * 		>
	 */
	// TODO: use same return type as in stupid implementation
	std::pair< bool, size_t > runFordBellman( size_t srcVid, size_t destVid, Cost * dist, GrEdgePtr * edgeFromParent )
	{
		size_t NOT_A_VERTEX_ID = GRAPH_SUPPORT_NS::Graph::NOT_A_VERTEX_ID;
		int n = gr.getVertexCount();
		Cost infCost = EdgeAttributes::getINFCost();
		Cost zeroCost = 0;
		std::fill( dist, dist + n, infCost );
		std::fill( edgeFromParent, edgeFromParent + n, NULL );
		Capacity epsCapacity = EPS;
		dist[srcVid] = zeroCost;
		bool updated = false;
		int iterationCounter = 0;
		bool negativeCycleExists = false;
		size_t lastUpdatedDestVid = NOT_A_VERTEX_ID;
		do
		{
			updated = false;
			for ( int vid = 0; vid < n; vid++ )
			{
				EdgeIterator stopeit = gr.getEdgsEndId( vid );
				for ( EdgeIterator eit = gr.getEdgsBegId( vid ); eit != stopeit; ++eit )
				{
					GrEdgePtr eptr = * eit;
					EdgeAttributes & attributes = eptr->getAttributes();
					if ( attributes.getRemainCapacity() < epsCapacity )
					{
						continue;
					}

					Cost ndist = dist[vid] + attributes.getCost();
					// TODO: use real double inf
					ndist = dist[vid] == infCost ? dist[vid] : ndist;
					size_t dVid = eptr->endId();
					if ( dist[dVid] > ndist + EPS )
					{
						dist[dVid] = ndist;
						edgeFromParent[dVid] = eptr;
						updated = true;
						lastUpdatedDestVid = dVid;
					}
				}
			}
			iterationCounter++;
			if ( iterationCounter > n )
			{
				negativeCycleExists = true;
				break;
			}
		} while ( updated );
		assert( ! negativeCycleExists || lastUpdatedDestVid != NOT_A_VERTEX_ID );
		bool pathFound = dist[destVid] < EdgeAttributes::getINFCost();
		assert( ! pathFound || edgeFromParent[destVid] != NULL );
		size_t negativeCycleChild = negativeCycleExists ? lastUpdatedDestVid : NOT_A_VERTEX_ID;
		return std::pair< bool, size_t >( pathFound, negativeCycleChild );
	}

	bool runPotentialDejkstra( size_t srcVid, size_t destVid, Cost * dist, Cost * fi, bool * visited, GrEdgePtr * edgeFromParent )
	{
		// TODO:
		size_t NOT_A_VERTEX_ID = GRAPH_SUPPORT_NS::Graph::NOT_A_VERTEX_ID;
		// TODO: move this typedef upper
		int n = gr.getVertexCount();
		std::fill( visited, visited + n, false );
		Cost infCost = EdgeAttributes::getINFCost();
		// TODO: eps - s in this class
		Cost epsCost = 1e-5;
		Cost mepsCost = -1e-5;
		Cost zeroCost = 0;
		std::fill( dist, dist + n, infCost );
		dist[srcVid] = 0;
		std::fill( edgeFromParent, edgeFromParent + n, NULL );
		Capacity epsCapacity = EPS;
		dist[srcVid] = zeroCost;
		int iterationCounter = 0;
		for ( int i = 0; i < n; i++ )
		{
			int minVid = NOT_A_VERTEX_ID;
			Cost minCost = infCost;
			for ( int i = 0; i < n; i++ )
			{
				if ( minCost > dist[i] && ! visited[i] )
				{
					minCost = dist[i];
					minVid = i;
				}
			}
			if ( minVid == NOT_A_VERTEX_ID )
			{
				break;
			}
			assert( ! visited[minVid] );
			visited[minVid] = true;
			// TODO: common stype of stopeit;
			EdgeIterator stopeit = gr.getEdgsEndId( minVid );
			for ( EdgeIterator eit = gr.getEdgsBegId( minVid ); eit != stopeit; ++eit )
			{
				GrEdgePtr eptr = * eit;
				Capacity remainCapacity = eptr->getAttributes().getRemainCapacity();
				if ( remainCapacity < epsCapacity )
				{
					continue;
				}
				size_t dvid = eptr->endId();
				Cost additCost = eptr->getAttributes().getCost() + fi[minVid] - fi[dvid];

				assert( additCost >= mepsCost );
				Cost ncost = minCost + additCost;
				if ( dist[dvid] > ncost + epsCost )
				{
					dist[dvid] = ncost;
					edgeFromParent[dvid] = eptr;
				}
			}
			iterationCounter++;
			assert( iterationCounter <= n );
		}
		return visited[destVid];
	}

	void pushViaPath( size_t startVid, size_t finishVid, GrEdgePtr * edgeFromParent, Capacity desiredCapacity )
	{
#ifndef NDEBUG
		std::vector< char > visitedBuf( gr.getVertexCount(), ( char ) false );
		bool * visited = ( bool * ) & visitedBuf[0];
#endif
		size_t vid = finishVid;
		Capacity minCapacity = desiredCapacity;
		do
		{
			assert( ! visited[vid] ); assert( visited[vid] = true );
			GrEdgePtr & eptr = edgeFromParent[vid];
			assert( eptr != NULL );
			Capacity remainCapacity = eptr->getAttributes().getRemainCapacity();
			if ( minCapacity > remainCapacity )
			{
				minCapacity = remainCapacity;
			}
			vid = eptr->begId();
		} while ( vid != startVid );

		vid = finishVid;
		Cost fee = 0;
		do
		{
			GrEdgePtr & eptr = edgeFromParent[vid];
			assert( eptr != NULL );
			EdgeAttributes & attrs = eptr->getAttributes();
			attrs.addFlow( minCapacity );
			fee += attrs.getCost() * minCapacity;
			vid = eptr->begId();
		} while ( vid != startVid );

		Capacity zeroCapacity = 0;
		totalFlow += startVid == finishVid ? zeroCapacity : minCapacity;
		totalCost += fee;
	}

	size_t discoverNegativeCyclePart( size_t negativeCycleChild, GrEdgePtr * edgeFromParent, int * visited, int & curVisitedTrue )
	{
		curVisitedTrue++;
		size_t vid = negativeCycleChild;
		while ( visited[vid] != curVisitedTrue ) {
			visited[vid] = curVisitedTrue;
			GrEdgePtr & eptr = edgeFromParent[vid];
			assert( eptr != NULL );
			vid = eptr->begId();
		}
		assert( visited[vid] == curVisitedTrue );
		return vid;
	}

private:
	static const double EPS = 1e-5;

private:
	Graph & gr;
	Capacity totalFlow;
	Cost totalCost;
};

}
}
}

#endif /* MINCOSTFLOW_H_ */
